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Synopsis Thermal limits in ectotherms may arise through a mismatch between supply and demand of oxygen. At higher 
temperatures, the ability of their cardiac and ventilatory activities to supply oxygen becomes insufficient to meet their 
elevated oxygen demand. Consequently, higher levels of oxygen in the environment are predicted to enhance tolerance of 
heat, whereas reductions in oxygen are expected to reduce thermal limits. Here, we extend previous research on thermal 
limits and oxygen limitation in aquatic insect larvae and directly test the hypothesis of increased anaerobic metabolism 
and lower energy status at thermal extremes. We quantified metabolite profiles in stonefly nymphs under varying tem- 
peratures and oxygen levels. Under normoxia, the concept of oxygen limitation applies to the insects studied. Shifts in the 
metabolome of heat-stressed stonefly nymphs clearly indicate the onset of anaerobic metabolism (e.g., accumulation of 
lactate, acetate, and alanine), a perturbation of the tricarboxylic acid cycle (e.g., accumulation of succinate and malate), 
and a decrease in energy status (e.g., ATP), with corresponding decreases in their ability to survive heat stress. These 
shifts were more pronounced under hypoxic conditions, and negated by hyperoxia, which also improved heat tolerance. 
Perturbations of metabolic pathways in response to either heat stress or hypoxia were found to be somewhat similar but 
not identical. Under hypoxia, energy status was greatly compromised at thermal extremes, but energy shortage and 
anaerobic metabolism could not be conclusively identified as the sole cause underlying thermal limits under hyperoxia. 
Metabolomics proved useful for suggesting a range of possible mechanisms to explore in future investigations, such as the 
involvement of leaking membranes or free radicals. In doing so, metabolomics provided a more complete picture of 
changes in metabolism under hypoxia and heat stress. 



Introduction 

Thermal limits in ectotherms may arise through a 
mismatch between supply and demand of oxygen 
(Winterstein 1905; Portner 2001). As temperatures 
increase, there is an increase in both the environ- 
mental availability of oxygen and organismal 
demand for oxygen, with the latter tending to 
increase more rapidly (Verberk et al. 2011). As a 
result, the ratio of oxygen supply relative to 
demand decreases and shortages of oxygen may 
arise in warmer waters. Consequently, cardiac and 



ventilatory activities of ectotherms may be insuffi- 
cient to meet their elevated demand for oxygen at 
the higher temperatures (Frederich and Portner 
2000). The resulting whole-animal drop in aerobic 
scope induces a shift from aerobic to anaerobic me- 
tabolism to maintain energy status (Portner 2006). 

Evidence for temperature-dependent oxygen limi- 
tation has not been forthcoming for terrestrial insects 
(Klok et al. 2004; Stevens et al. 2010); most of the 
evidence for oxygen limitation at thermal extremes 
to date comes from a variety of marine taxa 
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(Portner 2001). In tracheated arthropods, oxygen lim- 
itation may be unlikely due to the efficiency and plas- 
ticity of tracheal systems in supplying oxygen directly 
to metabolically active tissues. Recently, Verberk and 
Bilton (2011) found that better oxygenation of the 
medium improves heat tolerance in the aquatic 
nymphs of the stonefly Dinocras cephalotes (Curtis 
1827). These nymphs have closed tracheal systems 
and rely on diffusion across their tegument for gas 
exchange and are hence limited in their regulatory 
capacity. Although these findings support the idea 
that a mismatch between external oxygen supply and 
internal oxygen demand can set thermal limits in aqua- 
tic insects they do not conclusively show that thermal 
limits are caused by a loss of aerobic scope of the whole 
animal. This requires demonstrating the onset of an- 
aerobic metabolic pathways at thermal extremes and a 
modulating effect of oxygen on the onset of these 
pathways. Here, we extend previous research on ther- 
mal limits and oxygen limitation in aquatic insect 
larvae and investigate oxygen limitation directly by 
looking at the onset of anaerobic metabolism. 

Generally, anaerobic metabolism is not intensively 
studied in insects (Miiller et al. 2012) and specific 
anaerobic end products are unknown for many in- 
sects, including this species. Therefore, rather than 
measuring a small set of metabolites that are typically 
associated with anaerobic metabolism (e.g., lactate, 
succinate, and alanine) but risk missing the actual 
end products, we chose to use an untargeted metabo- 
lomic approach. Metabolomics is the now well-estab- 
lished technology concerned with the study of 
naturally occurring, low-molecular-weight organic 
metabolites within a biological specimen. Environ- 
mental metabolomics can be defined simply as the 
application of this technology to characterize the in- 
teractions of organisms with their environment 
(Malmendal et al. 2006; Bundy et al. 2009; Colinet 
et al. 2012; Viant and Sommer 2013). This approach 
has considerable potential for characterizing the 
responses of aquatic organisms to natural and anthro- 
pogenic stressors (Viant 2007). Here, we used 
metabolomics in stonefly nymphs for the first time 
to discover (1) which metabolites are accumulating 
and if they indicate the onset of anaerobic metabolism 
and (2) whether the pattern in metabolites supports 
the hypothesis of increased anaerobic metabolism and 
lower energy status at thermal extremes. 

Methods 

Maintenance of animals 

The life-cycle of the stonefly D. cephalotes lasts 3 
years, and most of this time is spent as aquatic 



nymphs. Aquatic perlid nymphs such as D. cephalotes 
are among the largest European freshwater carnivo- 
rous invertebrates, feeding primarily at night 
(Malmqvist and Sjostrdm 1980). Nymphs pass 
through 15-21 instars, and growth is slower for 
males than females, resulting in pronounced sexual 
dimorphism (Frutiger 1987). Maximum growth of 
nymphs occurs in spring-summer. Adult flight and 
oviposition occur chiefly in May and June. 

Aquatic nymphs were collected from the River 
Dart, Devon, UK at the same site in (early) spring 
in 2010 and (late) spring in 2011. Higher air tem- 
peratures and less rainfall were recorded for SW 
England in 2011 compared with 2010 for the 
period prior to animal collection (January- April). 
Nymphs were housed individually in separate cham- 
bers that were placed in a flow-through aquarium 
(lOlmin -1 ), fed with artificial pond water (ASTM 
1980), buffered, and diluted to reflect the pH and 
conductivity of the field site (pH 6.4-6.6, 70- 
150 ).iS cm" 1 ). Nymphs were fed chironomid larvae 
and maintained in the laboratory at 10±1°C in a 
12 L:12 D regime. We could maintain nymphs under 
these conditions for extended periods of time 
(months), and all nymphs were acclimated for at 
least 7 days to laboratory conditions before trials 
were undertaken to measure endpoints used to 
assess critical temperatures. Nymphs could not be 
sexed reliably, and therefore, we could not account 
for differences between sexes, which may have in- 
creased the variance in our results. 

Oxygen consumption 

Oxygen consumption was measured for each larvae 
at 5, 10, and 15°C using closed glass respiration 
chambers of 67.5-68.5 ml. The chambers were im- 
mersed in a temperature controlled bath (±0.1°C) 
and stirred using submersible magnetic stirrers to 
ensure mixing of water. Respiration chambers were 
fitted with a fine nylon mesh forming a false bottom 
to prevent contact between the larva and the mag- 
netic stirrer bar. Individuals were allowed to accli- 
mate for 60 min before the chambers were closed and 
left for at least 60 min (60-140 min). Oxygen content 
was measured before closing the chambers and after 
the incubation period using an oxygen electrode 
(1302 Oxygen Electrode, Strathkelvin Instruments) 
connected to a calibrated meter (Oxygen Meter 
929, Strathkelvin Instruments). On average, larvae 
depleted oxygen levels to 94% of the initial values 
(minimum 82%), and oxygen consumption was ex- 
pressed as ug 0 2 tT 1 g fresh weight -1 . 
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Critical temperatures 

To assess critical temperatures, we employed the 
same methods as previously described (Verberk and 
Bilton 2011; Verberk and Calosi 2012). Briefly, indi- 
vidual nymphs were placed in five parallel flow- 
through chambers (70 x 70 x 30 mm; flow rate 
0.016 Is -1 ), and water was supplied to these cham- 
bers by gravity from a 25-1 header tank after having 
passed through a tubular counter-current heat 
exchanger. Water in the header tank was of the 
same composition as that used to maintain animals 
and was bubbled with a mixture of 20% 0 2 and 80% 
N2, obtained using a gas-mixing pump (Wosthoff, 
Bochum, Germany). Individuals were left resting 
for 1 h at the equilibration temperature of 10°C, 
after which temperature in the experimental cham- 
bers was increased to 0.25°Cmin _ , using a Grant R5 
water bath with a GP200 pump unit (Grant 
Instrument Ltd, Cambridge, UK), connected to the 
heat exchanger. Temperatures were logged using a 
HH806AU digital thermometer (Omega 
Engineering Inc., Stamford, CT). 

Critical maximal temperatures (CT max ) were re- 
corded as the point at which animals no longer 
showed any body movement or muscular spasms. 
Animals that were transferred at this point to fully 
oxygenated water of 10°C recovered (i.e., they sur- 
vived and performed equally well in subsequent 
trials). At lower temperatures, larvae initiated in 
repeated swimming behavior (at about 29°C; inter- 
preted as attempts to escape experimental condi- 
tions) and fell upon their backs (at about 30°C), 
and an onset of spasms was observed (at about 
31°C). CT max was assessed at normoxia, hypoxia, 
and hyperoxia. Different levels of oxygenation were 
achieved by changing the 0 2 -N 2 gas mixture 
obtained using the gas-mixing pump (Wosthoff, 
Bochum, Germany). The gas mixture was adjusted 
10 min after placing the animals in the small flow- 
through chambers, to allow for gradual exposure to 
hypoxic and hyperoxic conditions during the 1 h 
resting period. During this period, oxygen levels in 
the outflow water from the chambers were measured 
repeatedly, to verify that the oxygen levels had stabi- 
lized. Trials to assess CT max were performed in 2010 
and 2011. In 2010, mild hypoxic and hyperoxic con- 
ditions were used (14 and 36kPa, respectively). In 
2011, more extreme hypoxic and hyperoxic condi- 
tions were used (5 and 60kPa, respectively). 

Design of the study of metabolomics 

At the end of the trial, when animals reached the 
critical temperature, they were freeze-clamped in 



liquid nitrogen. These samples constituted three 
treatments, one for each level of oxygenation (hyp- 
oxia, normoxia, and hyperoxia). Two additional 
treatments comprised (1) animals freeze-clamped in 
liquid nitrogen before commencement of a heating 
trial (i.e., control animals at 10°C) and (2) animals 
subjected to a warming trial under hyperoxic condi- 
tions, but freeze-clamped in liquid nitrogen after 
reaching the critical temperature observed under 
hypoxic conditions (i.e., before reaching their critical 
temperature). All animals thus frozen were stored at 
— 85°C until further analysis. So in total there are 
five treatments: 

(1) Hypoxia at CT max . 

(2) Normoxia at CT max . 

(3) Hyperoxia at CT max . 

(4) Hyperoxia at CT hypoxia (treatment 1)- 

(5) Control (acclimation at 10°C, 20kPa). 

The first three treatments were all at the thermal 
limit and are therefore expected to be at a tempera- 
ture at which oxygen limitation has set in. In this 
sense, they could be said to have been standardized 
for this effect, even though the absolute temperatures 
varied in accordance with the level of oxygenation. 
In the last two treatments, animals were freeze- 
clamped in liquid nitrogen before reaching their 
thermal limit, with the control animals experiencing 
no heat stress, whereas the animals under hyperoxia 
in treatment 4 did experience the same heat stress as 
the first treatment, but (presumably) before oxygen 
limitation set in. These different conditions across 
these treatments allow three analyses that provide 
insight in the metabolites that are indicative of 
heat stress, oxygen limitation, or both. The first 
comparison is between nymphs at their critical tem- 
perature (treatments 1-3) and control animals (treat- 
ment 5). The metabolites differentiating these two 
groups provide a test of the hypothesis that thermal 
tolerance is limited by a shortage of oxygen and 
energy (as indicated by anaerobic end products and 
low ATP levels). The second comparison is between 
nymphs at their critical temperature but at different 
levels of oxygen (treatments 1-3). This allows a test 
of the hypothesis that metabolic profiles of animals 
at their critical temperature are comparable irrespec- 
tive of the oxygen conditions as they should all be 
experiencing oxygen limitation. The third compari- 
son between nymphs under hypoxia and hyperoxia, 
controlling for temperature (treatments 1 and 4), 
provides a test as to whether the increased heat tol- 
erance under hyperoxia is associated with an allevi- 
ation or postponement of anaerobic metabolism. 
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Extraction of metabolites 

Metabolites were extracted from each biological 
sample using a two-step methanol/water/chloroform 
protocol (Wu 2008). In brief, per milligram of tissue, 
1 pi methanol (high-performance liquid chromatog- 
raphy [HPLC] grade) and 0.2 pi water (HPLC grade) 
were added to each sample (one freeze- clamped 
stonefly nymph) prior to homogenization, using a 
Precellys-24 ceramic bead-based homogenizer 
(Stretton Scientific Ltd, UK). Thereafter, per milli- 
gram of tissue, 0.5 pi water (HPLC grade) and 1 pi 
chloroform (pesticide analysis grade) were added to 
each sample, which was vortexed and centrifuged 
(4000 rcf, lOmin), yielding an upper (polar) and 
lower (non-polar) fraction for each sample. Each 
polar fraction was collected into a 1.5-ml 
Eppendorf tube, and the non-polar fraction into a 
1.8-ml glass vial. Of the polar extract, 5 pi was taken 
for mass spectrometry and 50 pi for NMR spectros- 
copy. All samples were dried using a centrifugal con- 
centrator (Thermo Savant, Holbrook, NY) and 
stored at — 80°C prior to analysis. 

FT-ICR mass spectrometry and spectral processing 

Immediately prior to mass spectrometric analysis, 
polar fractions were resuspended in 40 pi methanol/ 
water (4:1) with 20 mM ammonium acetate for neg- 
ative ion analysis (the non-polar samples were not 
analyzed), and 5 pi each was taken to pool a quality 
control sample. For initial optimization of concen- 
tration, test samples were further diluted 1:1, 1:3, 1:7, 
and 1:15; following this, the main samples were not 
further diluted. Samples were centrifuged (lOmin, 
22,000 rcf, 4°C) to remove any particulate matter. 
Fourier transform ion cyclotron resonance (FT- 
ICR) direct infusion mass spectrometry was per- 
formed using an LTQ FT Ultra (Thermo Fisher 
Scientific, Bremen, Germany) equipped with a 
chip-based direct infusion nanoelectrospray ion 
source (Triversa, Advion Biosciences, Ithaca, NY). 
Each sample was analyzed in triplicate from a 96- 
well plate using the selected ion monitoring (SIM)- 
stitching method as adapted for the LTQ FT Ultra 
(Weber et al. 2011) from m/z 70 to 590 in both ion 
modes. Raw mass spectral data were processed using 
the SIM-stitching algorithm (Southam et al. 2007; 
Payne et al. 2009). Missing values were filled into 
the resulting data matrix using the KNN algorithm 
(Hrydziuszko and Viant 2012) in an in-house R 
script, whereas PQN normalization (Dieterle et al. 
2006) was performed using again an in-house 
Matlab script. This normalized matrix was trans- 
formed using the generalized logarithm (Parsons 



et al. 2007) to stabilize the variance across the thou- 
sands of measured peaks prior to analysis with mul- 
tivariate statistics (see below). 

NMR spectroscopy and spectral processing 

The dried polar extracts were resuspended in 650 pi 
sodium phosphate buffer in D 2 0 (0.1 M, pH 7.0) 
containing 0.5 mM sodium 3-trimethylsilyl-2,2,3,3- 
d4-propionate, which serves as an internal chemical 
shift standard. The samples were centrifuged, and 
then 600 pi of each supernatant was transferred to 
a 5 -mm NMR tube and analyzed immediately 
using a DRX-500 NMR spectrometer (Bruker 
Biospin, Coventry, UK), equipped with a cryoprobe 
and operated at 500.18MHz (at 300 K). One-dimen- 
sional (ID) : H NMR spectra were obtained as de- 
scribed by Hines et al. (2007). Spectra were obtained 
using excitation sculpting for water suppression 
(Hwang and Shaka 1995) and using an 8.4 ps (60°) 
pulse, 6 kHz spectral width, and a 2.5 s relaxation 
delay with water presaturation. A total of 64 tran- 
sients were collected into 16,348 data points, requir- 
ing a 4.5-min acquisition time. Datasets were zero- 
filled to 32,768 points, before line-broadenings of 
0.5 Hz were applied prior to Fourier transformation. 

To maximize metabolite discrimination two-di- 
mensional (2D) : H J-resolved (JRES) NMR spectra 
were also acquired (Viant 2003), being processed ac- 
cording to Hines et al. (2007). 2D JRES spectra were 
acquired for each sample using 16 transients per in- 
crement, for 16 increments, which were collected 
into 16,000 data points with spectral widths of 
6 kHz in F2 (chemical shift axis) and 50 kHz in Fl 
(spin-spin coupling constant axis). A 4.0-s relaxation 
delay was employed resulting in a total acquisition 
time of 24min. Datasets were zero-filled in Fl, the 
F2 dimension was then multiplied by a SEM window 
function using 0.5 Hz line broadening while the Fl 
dimension was multiplied by a sine-bell window 
function, all prior to Fourier transformation. 2D 
JRES spectra were tilted by 45°, symmetrized about 
Fl, and calibrated using TopSpin (Bruker Biospin). 
Data were exported as the ID skyline projections of 
2D JRES spectra and converted to a format for mul- 
tivariate statistical analyses using NMRLab (Giinther 
et al. 2000). 

Statistical analysis and identification of metabolites 

Principal components analysis (PCA) was used ini- 
tially to assess the overall metabolic differences be- 
tween the sample groups in an unbiased manner, 
using the PLS_Toolbox (version 5.5.1, Eigenvector 
Research, Manson, WA) within Matlab (version 
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7.8; The Maths Works, Natick, MA). Supervized mul- 
tivariate analyses were performed using partial least- 
squares discriminant and regression approaches 
(PLS-DA and PLS-R, respectively) using the PLS 
Toolbox, with internal cross-validation (Venetian 
blinds) and permutation testing (1000 tests each) 
using in-house Matlab scripts. Univariate statistical 
analyses were used to confirm the significance 
of changes in individual mass-spectral signals. 
Specifically, analysis of variance was conducted 
using in-house R scripts (with a false discovery rate 
of 5% to correct for multiple hypothesis testing) 
(Benjamini and Hochberg 1995). For FT-ICR mass- 
spectrometry data, identification of putative empiri- 
cal formulae was achieved using Mi-Pack software 
(Weber and Viant 2010), the KEGG database 
(http://www.genome.jp/kegg/download/), and calcu- 
lations in Xcalibur 2.0.7 (Thermo Fisher Scientific). 
For NMR data, peaks were identified and integrated 
using in-house peak-picking algorithms (J. Byrne, C. 
Ludwig, J. M. Easton, S. He, and M. R. Viant, in 
preparation). 

Results 

Ecophysiology of stonefly nymphs 

Oxygen consumption by stonefly nymphs increased 
with increasing temperatures and body size, with 
temperature effects being stronger in larger bodied 
animals (F lj207 = 27.16; P<0.0001; Fig. 1). Overall, 
oxygen consumption was comparable for nymphs 
collected in 2010 and 2011 (F lj207 = 1.51; P=0.221; 
Supplementary Table SI), although slightly stronger 
responses to temperature were found for nymphs 
collected in 2011 (F h207 = 5.53 7; P= 0.0196; 
Supplementary Table SI). The critical thermal 
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Fig. 1 Oxygen consumption of stonefly nymphs at 5, 10, and 
15°C shown separately for small (<260mg) and large nymphs 
(>260mg). Values in brackets refer to the numbers of measure- 
ments. Note that mass specific rates of oxygen consumption are 
shown; absolute rates of oxygen consumption were higher for 
larger animals. 



maxima (CT max ) observed by behavioral responses 
was consistently influenced by the level of oxygen 
(•Fi,77 = 31.22; P<0.0001), suggesting oxygen limita- 
tion at thermal extremes (Fig. 2 and Supplementary 
Table S2). Interannual differences in CT max were 
found (Pi )77 = 9.88; P= 0.0024), with nymphs col- 
lected in 2011 displaying lower CT max at normoxia 
(36.1°C in 2010 versus 35.0°C in 2011). In addition, 
there was a modest effect of body size, with larger 
sized individuals displaying lower heat tolerance 
(£=-3.53, t lj76 = -2.052, P=0.044). Nevertheless, 
treatment effects of the level of oxygen predomi- 
nated; relative to normoxia, hypoxia decreased 
CT max by 2.9°C (at 14kPa) and 7.8°C (at 5kPa), 
whereas hyperoxia increased CT max by 1.5°C (at 
36kPa) and 1.9°C (at 60kPa). 

Metabolomics of stonefly nymphs 

Peaks from the FT-ICR mass spectra of homogenates 
of whole organisms were annotated with empirical 
formula(e) and putative metabolite names using the 
Mi-Pack software. Of the >3000 peaks detected, 
99.6% were assigned to one or more empirical for- 
mula^), representing the most extensive putative an- 
notation of the stonefly metabolome to date. 
Transformation mapping in conjunction with the 
KEGG database was then used to attempt to assign 
putative metabolite names to these empirical formu- 
lae (Weber and Viant 2010). In total, 196 peaks 
could be assigned to one or more compound 
names listed in the KEGG database. The metabolo- 
mics data presented here are based primarily on 
those derived from the FT-ICR mass spectra, and 
where NMR data were available for the same metab- 
olites, the measurements were compared across the 
two analytical approaches. Both methods yielded 
highly similar results (Supplementary Fig. SI), 
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Fig. 2 Critical thermal maxima (CT max ) at different levels of 
oxygen for stonefly nymphs collected in 2010 and 2011. Values in 
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which is interesting from a methodological point of 
view. In addition, the NMR data provided greater 
confidence in the identification of the metabolites, 
in particular when a signal in the mass spectra 
could originate from multiple putative metabolite 
names having the same empirical formula. For 
example, identification of lactate, succinate, fuma- 
rate, tyrosine, and alanine (Supplementary Table 
S3) is based on the fact that these metabolites were 
also identified in the NMR spectra. 

Preliminary data analysis using PCA showed that 
differences in metabolic profiles between the nymphs 
collected in 2010 and 2011 were subordinate to the 
differences observed across the different treatment 
(Supplementary Fig. S2). Therefore, the remainder 
of the analyses focuses on comparing the various 
treatments to elucidate metabolic responses to tem- 
perature and oxygen. 

Comparing the metabolites of nymphs at 10°C with 
those at their thermal limit 

To determine whether our metabolomics approach 
could distinguish between the metabolome of 
nymphs at their critical temperature (treatments 
1—3) and that of control animals (treatment 5), we 
conducted supervized multivariate classification 
(PLS-DA). The optimal PLS model comprised five 
latent variables (LVs), based on the minimization 
of the cross-validated classification errors, and 
accounted for 45.6% of the total variance in the 
metabolic dataset. The metabolic profile of stonefly 
nymphs reaching the critical thermal maxima 
(CT max ; >33°C), and control animals kept at 10°C 
were clearly distinguishable (Fig. 3). The cross- 
validated classification error was 15.9% with a signif- 
icance of P< 0.001 when compared with 1000 
random permutations. Furthermore, the metabolites 
important in separating these two groups did, to 
some extent, also differentiate within nymphs at 
their thermal limit according to oxygen level; 
nymphs under hyperoxia tended to be more similar 
in metabolic profile to control animals, and nymphs 
under hypoxia were least similar to control animals 
(Supplementary Table S3). The metabolites differen- 
tiating between both groups are associated with the 
generation of energy, involving sugar metabolism 
and anaerobic metabolism (e.g., lactate and succi- 
nate; see Fig. 4 and Supplementary Table S3). 
In addition, ATP levels were also lower in 
animals at their thermal limit. Uridine levels were 
elevated in tandem with lactate (Spearman rank cor- 
relation across all individuals: r— 0.734, P< 0.001). 
Furthermore, ATP levels were strongly related to 
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levels of L-arginine phosphate (Spearman rank cor- 
relation across all individuals: r= 0.405, P<0.001). 

Comparing the metabolism of nymphs reaching their 
thermal limit at different levels of oxygen 

To assess whether there were consistent differences in 
the metabolite profile of animals reaching their crit- 
ical thermal maxima but at different oxygen levels 
(treatments 1-3), we conducted supervized multivar- 
iate regression (PLS-R). Perturbations in the meta- 
bolic pathways were highly predictive of the oxygen 
conditions that nymphs experienced during the 
experiments. The model with four LVs accounted 
for 43.3% of the variance in the total metabolic data- 
set (x-block). Oxygen values (y) predicted from the 
metabolic profiles accurately matched those used in 
the experiment (x) (Fig. 5; partial least squares 
regression; y=0.981x; _R 2 = 0.944, Q 2 [cross-vali- 
dated] = 0.605, P<0.001 [1000 random permuta- 
tions]). This suggests that the metabolic responses 
to heat stress in the various oxygen treatments 
were distinct. The metabolites differentiating between 
oxygen conditions were again associated with the 
generation of energy, with anaerobic metabolism 
being more pronounced under hypoxia (treatment 1) 
(e.g., lactate and succinate) (see Fig. 4 and 
Supplementary Table S3). In addition, metabolism 
of amino acids was perturbed (e.g., alanine and 
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Fig. 5 Plot of measured versus PLS-R-predicted oxygen levels, 
showing a high correlation between observed and predicted 
oxygen levels (partial least squares regression; y = 0.981 x; 
R 2 = 0.944, Q 2 [cross-validated] = 0.605, P< 0.001 [1000 random 
permutations]). Plot is based on nymphs at their critical thermal 
maxima only (treatments 1-3). 



L-arginine phosphate) (see Fig. 4 and Supplementary 
Table S3). 
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Comparing the metabolism of nymphs under hypoxia 
and hyperoxia; controlling for temperature 

To test whether the observed improvement of heat 
tolerance under hyperoxic conditions (Fig. 2) is 
associated with alleviation or postponement of 
anaerobic metabolism, we compared the metabolic 
profile of nymphs reaching their thermal maximum 
under hypoxic conditions (treatment 1) with that of 
nymphs at the same temperature but under hyper- 
oxic conditions (treatment 4). The metabolic profiles 
were clearly different for both treatments as 
demonstrated by PLS-DA (Fig. 6). The optimal PLS 
model comprised three LVs, based on the minimiza- 
tion of cross-validated classification errors, and 
accounted for 36.7% of the total variance in the 
metabolic dataset (x-block). The cross-validated 
classification error was 14.4% with a significance of 
P< 0.001 when compared with 1000 random permu- 
tations. So there were consistent differences between 
hypoxia and hyperoxia in the perturbations in 
metabolic pathways in response to heat stress. 
Again, the metabolites differentiating between the 
two groups are associated with the generation of 
energy, involving both anaerobic metabolism (e.g., 
lactate and alanine) (see Fig. 4 and Supplementary 
Table S3) and amino-acid metabolism (e.g., 
tyrosine). 
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Fig. 6 PLS-DA scores plot, showing a clear separation in meta- 
bolic profile between nymphs under hypoxia (treatment 1) and 
hyperoxia (treatment 4) irrespective of the temperature (33.1 °C 
at 14kPa and 27.3 C at 5 kPa) (cross-validated class error rate of 
14.3%, P<0.001 [1000 random permutations]). The best model 
had three LVs in total, but for clarity only the first two axes are 
shown. The percentage of data variance explained by each latent 
variable is indicated. Ellipses are drawn by eye to improve 
visibility. 

Discussion 

The respiratory system of terrestrial arthropods con- 
sists of an open tracheal system, a highly branched 
system of tubes extending throughout the body and 
supplying oxygen directly to metabolically active tis- 
sues. The air-filled trachea opens to the external 
atmosphere and can be ventilated, which greatly 
facilitates gas exchange (Dejours 1981; Socha et al. 
2008). Consequently, insects with open trachea can 
regulate oxygen intake and maintain aerobic scope 
while avoiding oxygen toxicity under widely varying 
conditions (Hetz and Bradley 2005; Harrison et al. 
2006). These considerations could explain why 
evidence for oxygen limitation at thermal maxima 
has not been forthcoming in several tracheated 
arthropods (Klok et al. 2004; Stevens et al. 2010). 
In this study, we studied nymphs with closed tra- 
cheal systems and which rely on diffusion across 
their tegument for gas exchange. We extend previous 
research showing that results for terrestrial insects 
cannot be generalized to all insect taxa and life 
stages (Verberk and Bilton 2011). 

In the stonefly studied here, hyperoxia improved 
tolerance to heat whereas hypoxia reduced heat tol- 
erance (Fig. 1). More extreme levels of either hypoxia 
or hyperoxia did correspondingly change heat toler- 
ance to a greater extent, although the effects were not 



symmetrical: Hypoxia had larger effects than hyper- 
oxia (see also Verberk and Calosi 2012). In addition, 
heat tolerance differed among years by 1.1 °C under 
normoxia, showing that upper thermal limits are 
somewhat plastic in these aquatic nymphs, perhaps 
to a greater degree than in terrestrial insects 
(Overgaard et al. 2011; Hoffmann et al. 2013). 
Nymphs were collected from the same field site, 
and field conditions did not differ dramatically be- 
tween years (in fact conditions prior to collection 
were warmer and drier in 2011, contrasting with 
the lower upper thermal limits observed), suggesting 
that a greater plasticity in upper thermal limits might 
well be observed under more challenging conditions. 

Importantly, we directly investigated the link be- 
tween thermal limits and the expected onset of an- 
aerobic metabolism when oxygen limitation sets in. 
Higher levels of lactate and succinate clearly indicate 
the onset of anaerobic metabolism at thermal 
maxima in stonefly nymphs under normoxic condi- 
tions (Fig. 4). Higher levels of alanine, produced to- 
gether with lactate from pyruvate under conditions 
of limited oxygen, further support the onset of 
anaerobic metabolism. The pyruvate generated from 
glucose during glycolysis fits with the increased 
sugar-based metabolism we observed (e.g., pentitol 
phosphates). To reduce the accumulation of 
pyruvate, it can be transformed into acetylpyruvate 
using acetate (in addition to the transformation into 
lactate and alanine mentioned above). Although 
these specific end products were not known for 
D. cephalotes, they corroborate previous research on 
responses to hypoxia in Drosophila (Feala et al. 2007) 
and correspond to known by-products of anaerobic 
metabolism in other terrestrial insects (Hoback and 
Stanley 2001). Insects typically rely on the glycerol-3- 
phosphate shuttle to regenerate NAD+, which is im- 
portant for maintaining glycolysis (Gilmour 1961), 
thus explaining the higher levels of glycerol and glyc- 
erol-3-phosphate. Our findings in Dinocras show 
good correspondence to those by Malmendaal et al. 
(2006) who also include alanine, acetate, and glucose 
as important markers for heat stress in Drosophila. 

The higher dependency on glycolysis under 
oxygen-limited conditions is also evident from an 
overall disruption of the tricarboxylic acid (TCA) 
cycle with strong increases in oxaloacetate, malate, 
and fumarate, whereas levels of isocitrate decreased. 
Interestingly, the TCA cycle intermediate 2-oxogluta- 
rate (alfa-ketoglutaric acid) was unaffected by the 
treatments, but was strongly related to the body 
size of individuals, with larger animals having 
higher levels (£ = 0.843, f lj72 = 6.49, P< 0.001; 
Supplementary Fig. S3). Michaud et al. (2008) 
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studied in heat stress in the Antarctic midge and 
suggested that increased 2-oxoglutarate levels were 
indicative of a perturbed TCA cycle. In addition, 
2-oxoglutarate is involved in the formation of ala- 
nine from pyruvate in various metazoans, including 
sipunculid and polychaete worms (Miiller et al. 
2012). Larger animals also consistently showed 
reduced levels of acetylpyruvate (/3 = — 0.365, 
^,69 = 2.88, P= 0.0054). This suggests that metabolic 
pathways are size dependent, with the TCA cycle 
being perhaps more important relative to glycolysis 
in larger animals. Such size-dependent differences in 
metabolic pathways may underlie the size-dependent 
differences in both heat tolerance (Supplementary 
Table S2) and thermal sensitivity of metabolism 
found (Fig. 1; see also Verberk and Bilton 2011). 

The metabolic shifts toward anaerobic metabolism 
described above in response to oxygen limitation 
were not sufficient to compensate for the decrease 
in energy production via aerobic metabolism. ATP 
levels indicate a lower energy status in nymphs at 
their thermal maxima (treatments 1-3). In addition, 
Sacktor and Cochran (1957) found that in addition 
to hydrolysis of ATP to ADP, the hydrolysis of UTP 
to release the stored chemical energy involved split- 
ting all three phosphates, resulting in free uridine. 
Such generation of free uridine under conditions of 
ATP shortage could explain both the high levels of 
uridine under conditions at which energy status is 
compromised and the tight association between uri- 
dine and lactate accumulation. In addition, compar- 
ing nymphs at their thermal maxima (treatment 
1-3), arginine phosphate levels declined in parallel 
with the shift toward anaerobic metabolism. 
Arginine phosphate is a well-established marker of 
cellular energy status, donating a phosphate to 
ADP when the ATP pool is depleted (Uda et al. 
2006), explaining the positive correlation observed 
between ATP and arginine phosphate. 

These results clearly show that under normoxia, 
the concept of oxygen limitation applies to the 
insects studied. The onset of anaerobic metabolism 
and the decrease in energy status, which eventually 
limit survival of heat stress, can all be traced back to 
shifts in the metabolome of the stonefly nymphs. 
Comparing only nymphs at their thermal limit 
(treatments 1-3), we found an attenuation of these 
shifts in metabolites under hyperoxia and a potenti- 
ation under hypoxia. As nymphs reached their crit- 
ical thermal maxima at different temperatures under 
different oxygen conditions (Fig. 2), this indicates 
that perturbations of metabolic pathways in response 
to heat stress and hypoxia may be similar, but not 
identical (in which case changes in metabolism due 



to heat or hypoxia would be interchangeable and 
organisms at their thermal limit would have a very 
comparable metabolite profile). Possibly animals 
warmed under hypoxic conditions started to 
go into anaerobic metabolism quite some time 
before reaching their critical thermal maximum. 
Furthermore, when comparing the nymphs under 
hypoxia (treatment 1) and nymphs at the same tem- 
perature but under hyperoxic conditions (treatment 
4), it is clear that hyperoxia rescued these animals 
and negated the severe energy shortage and anaero- 
bic metabolism observed under hypoxia. Thus, the 
hypothesis that heat tolerance is limited by insuffi- 
cient energy as a consequence of oxygen limitation is 
strongly supported at normoxic and especially hyp- 
oxic conditions. However, for nymphs at their ther- 
mal limit under hyperoxia, anaerobic metabolism 
and associated changes were not very distinct. 
Shortage of energy resulting from oxygen limitation 
could still be a factor involved in limiting heat tol- 
erance, as there was a trend to lower ATP levels, 
higher alanine levels, and a general disruption of 
the TCA cycle (trend to higher fumarate and 
malate levels and lower isocitrate and oxaloacetate 
levels). Yet, the results allow an exploration to see 
whether other factors are involved in setting the level 
of heat tolerance. 

Alternative mechanisms underpinning heat toler- 
ance are related to membrane properties and protein 
denaturation (Feder and Hoffmann 1999; Portner 
2001). The consequences of heat on membranes are 
increased fluidity and membrane leakage, which will 
promote proton leaks and lead to an uncoupling of 
mitochondria. This will lower ATP production. As 
no strong decline in ATP was observed under hyper- 
oxia, even though this conformed to the highest 
temperature experienced by nymphs across all treat- 
ments, membrane fluidity is probably not the main 
driver. Nevertheless, we did find some evidence for 
disruption of membranes. Phosphocholine is an in- 
termediate in the synthesis of phosphatidylcholine in 
tissues, a major component of biological membranes. 
Phosphocholine levels were elevated in all animals at 
their thermal limit, suggesting some generic mem- 
brane response to heat stress. Subsequent experi- 
ments would be necessary to confirm the 
importance of reorganization of membranes. Dena- 
turation of proteins is another potentially important 
mechanism underlying heat tolerance (Somero 
1995). In our study, we did see disruption of 
amino acid metabolism and nucleotide metabolism, 
but this was most likely due to active catabolism of 
proteins to generate energy, rather than protein de- 
naturation. Thermal limits in terrestrial insects are 
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usually higher, making protein denaturation more 
likely (see also Verberk and Bilton 2011; Verberk 
and Calosi 2012). 

Specifically for what may cause heat coma under 
hyperoxia, we also have to consider the effects of 
radical oxygen species, which may be more prevalent 
during metabolic disruption. Levels of 2-hydroxybu- 
tyrate were elevated relative to control animals, espe- 
cially under hypoxia and normoxia (treatments 1 
and 2). Higher levels could indicate an upregulation 
of protection against the action radical oxygen spe- 
cies, as 2-hydroxybutyrate is released as a byproduct 
when cystathionine is cleaved to cysteine that is in- 
corporated into the antioxidant glutathione. This fits 
with the strong disruption of the TCA cycle under 
hypoxia and normoxia. However, we found no com- 
pelling evidence for damage by radical oxygen species 
that set levels of heat tolerance levels under hyper- 
oxic conditions. 

Finally, there is the possibility that neuronal motor 
control is impacted in the nymphs before the onset 
of anaerobic metabolism and the associated organ- 
ism-wide shortage of energy. Dawson-Scully et al. 
(2010) found that activation of the Protein kinase 
G (PKG) pathway increased survival in adult 
Drosophila, but decreased locomotory function. For 
aquatic damselfly nymphs, heat hardening was pre- 
viously shown to increase survival but decrease loco- 
motory control, possibly via activation of the PKG 
pathway (Verberk and Calosi 2012). Therefore, this 
pathway could uncouple the behavioral responses 
observed (Fig. 2) and the biochemical responses in 
the metabolome. 

In conclusion, this study provides broad support for 
the hypothesis of oxygen limitation; anaerobic metab- 
olites accumulated in nymphs at their thermal limit 
under hypoxic and normoxic conditions. The metabo- 
lomics approach showed that energy shortage and 
anaerobic metabolism could not be conclusively iden- 
tified as the cause underlying thermal limits under 
hyperoxia. Metabolomics proved useful for suggesting 
a range of possible mechanisms to explore in future 
investigations, such as the involvement of leaking 
membranes or free radicals. In doing so, metabolomics 
provided a more complete picture of changes in me- 
tabolism under hypoxia and heat stress. 
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